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ABSTRACT 

We have studied deuterium fractionation on interstellar grains with the use of an exact 
method known as the direct master equation approach. We consider conditions perti- 
nent to dense clouds at late times when the hydrogen is mostly in molecular form and 
a large portion of the gas-phase carbon has already been converted to carbon monox- 
ide. Hydrogen, oxygen, and deuterium atoms, as well as CO molecules, are allowed to 
accrete onto dust particles and react there to produce various stable molecules. The 
surface abundances, as well as the abundance ratios between deuterated and normal 
isotopomers, are compared with those calculated with the Monte Carlo approach. We 
find that the agreement between the Monte Carlo and the direct master equation 
methods can be made as close as desired. Compared with previous examples of the 
use of the direct master equation approach, our present method is much more efficient. 
It should now be possible to run large-scale gas-grain models in which the diffusive 
dust chemistry is handled "exactly" . 
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with the differential rate equations used to model the 
gas-phase che mistry. A comparison of the two approache s is 
to be found in lStantcheva. Shematovich fc Herbsd ((2002) . 

We have, so far, successfully applied the direct master 
equation approach to a system consisting of H and O atoms 
and CO molecules accreting onto a grain s urface, and react- 
ing to form a variety of reaction products jStantcheva et all 
2002). In this paper, we report an expansion of the system 
to include D atoms accreting onto the surface and reacting 
with the other reactive surface species to produce deuter- 
at ed isotopome r s. Th e system is identical to that studied 
in ICaselli et al.l i2002fl with modified rate equations and a 
Monte Carlo approach. A somewhat smaller n etwork of re- 
ac tions was utilised b y lcharnlev et alJ (ll997Tl . As opposed 
to lCaselli et"afl l|2002|h we concentrate on the doubly, triply 
and quadruply-deuterated species which we believe are likely 
candidates for obser vation in highly fr actionated low-mass 
protostellar sources dParise et alJl2002D . The main purpose 
of this paper; however, is to show how the direct master 
equation method can be used efficiently. 

The paper is organised as follows. First, in Section [5] a 
short description of the system under investigation is given. 
The theory, including the formulae used in the calculations, 
is given in Section [2] Section [I] presents the results, and is 
followed by a discussion in Section 



1 INTRODUCTION 

Rate equations have long been used in modelling the 
surface chemistry that occur s on cold dust p a rticles 
in the interstellar medium Pickles fc Williams! Il977t 
[ffesegw^^eAs^^^eun^E^S 1 It has been pointed out 
iCharnlev. Tielens fc Rod gers 1997]) that, in the case when 
the average number of reactive species on a grain is suf- 
ficiently small, the rate equations may not be suitable to 
describe the system, since they do not take into considera- 
tion its discrete n ature. One alternative is to use a M onte 
Carlo approach llTielens fc Hagenl Il982t ICharnlev et alJ 
ll997tlCaselli et al.ll2002l:ICharnlevll200lL 

Recently ano t her approac h to this problem 

JCreen et all l200lt iBiham et all l200ll) . has been pro- 
posed. Like the Monte Carlo method, this approach is 
based on a master equation. But, unlike the Monte Carlo 
technique, the master equation is directly converted into 
differential equations. In particular, the rate equations for 
the population of highly reactive species with small surface 
abundances are replaced by sets of differential equations 
for the probability that certain number of these atoms 
or molecules exist on the surface at the same time. One 
advantage of this new exact approach is its easy coupling 
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Table 1. H, O, and CO gas-phase abundances (cm -3 ) utilised 



Abundance n 


Low 


Intermediate 


High 


H 


1.15 


1.15 


1.10 


O 


0.09 


0.75 


7.0 


CO 


0.075 


0.75 


7.5 


D 


0.3 


0.3 


0.3 



2 THE CHEMICAL NETWORK AND 
PHYSICAL CONDITIONS 

We have considered a system in which H, O, and D atoms 
and CO molecules accrete onto the grain surface and re- 
act to form the stable molecules H2, HD, D2, H2O, HDO, 
D 2 0, H2CO, HDCO, D2CO, CH3OH, CH3OD, CH2DOH, 
CH2DOD, CHD2OH, CHD2OD, CD3OH, CD3OD, 2 , 
C0 2 , and the highly reactive radicals OH, OD, HCO, DCO, 
H3CO, H2DCO, HD2CO, and D3CO. The system pertains 
to old dense interstellar clouds, in which hydrogen is mostly 
molecular and carbon is mostly in the form of CO. The cal- 
culations were performed mainly at a temperature T=10 K, 
for a period of 10 4 years, and fixed gas-phase abundances 
of the accreting species, an approximation justified by the 
relatively short period of evolution of the surface chemistry 
compared with the time scale of the gas-phase chemical pro- 
cesses. At times even earlier than 10 4 yr, a steady-state con- 
dition is reached in which the surface populations of the 
stable species grow linearly with increasing time while the 
reactive species have fixed populations. The number of avail- 
able sites per grain, which helps to define the diffusion rates, 
was taken to be 10 6 ; this number refers to so-called classical 
grains with a size of 0.1 /im. 

Three different sets of values for the gas-phase densi- 
ties of H, O, CO, were used throughout the calculations. 
Given in Table Q these sets of values are referred to as 
the low, intermediate, and high density cases because they 
were derived from gas- phase models at total densities of 10 3 , 
10 4 , and 10 s cm -3 , respectively. In all cases, the concentra- 
tion of H is near 1 cm -3 . Unless stated to the contrary, all 
results are for a gas-phase abundance for D of 0.3 cm" 3 . 
This very high value relative to the atomic H abundance 
is presumably produced via fractionation in the cold gas, 
although cur rent gas-phase models cannot ea sily reproduce 
such a value (iRoberts. Herbst fc Millarll2002ft . Variations in 
the atomic deuterium abundance are also considered. 

In calculating the accretion, desorption, and diffusion 
ra tes of the surface speci es, we follow ed the methods use d 
in iHaseeawa et alJ <ll992t) - see also ICaselli et aD feOOSft : 
ICaselli. Haseeawa fc Herbst! <ll996t) - and used their values 
for the parameters necessary for the calculations. All par- 
ticles were considered to diffuse over the surface solely via 
thermal hopping except for H and D atoms, for which quan- 
tum tunnelling was also considered. The rates used for dif- 
fusion are the so-called fast rates iStantcheva et aD 120021) 
because these magnify the differences between exact and 
approximate methods for studying diffusive surface chem- 
istry. Table [5] shows the accretion, evaporation, and diffu- 
sion rates over an entire grain for the accreting species in 
the model. Note that at the temperatures considered, the 
other heavy species in the model do not diffuse or evaporate 
at non-negligible rates. 



Table 2. Assorted rates for selected species at 10 K 



Species 


^acc 


(cm 3 s 4 ) 


t" 1 (s- 1 ) t" 1 

■•cvap V 13 j cliff 




H 


1.45(-5) 


1.88(-3) 5.14(4-4) 


D 


1.02(-5) 


1.67(-4) 3.92(4-2) 


O 


3.62(-6) 


2.03(-23) 4.24(-5) 


CO 


2.73(-6) 






Table 3. Surface reactions 


in the H,0,CO model. 


Number 




Reaction 


Ea (K) a 


1 


H 4 


H 


rj 

> H2 




2 


H 4 


O 


— > OH 




3 


H 4 


- OH 


TT / \ 

> H2U 




4 


H 4 


- CO 


— > HCO 


2000 


5 


H 4 


- HCO 


— ► H2CO 




6 


H 4 


H 2 CO 


— ► H3CO 


2000 


7 


H 4 


- H3CO 


/-ITT TTTT 




8 


H 4 


- D 


— > HD 




9 


H 4 


OD 


— ► HDO 




10 


H 4 


DCO 


— > HDCO 




ii 


H 4 


- HDCO 


— > H 2 DCO 




12 


H 4 


- D 2 CO 


— > HD2CO 


1925 


13 


H 4 


- H 2 DCO 


— > CH 2 DOH 




14 


H 4 


HD2CO 


— > CHD2OH 




15 


H 4 


- D3CO 


— > CD3OH 




1(5 


O 4 


- O 


— > Oa 




17 


O 4 


- CO 


— > CO2 


1000 


18 


O 4 


- HCO 


— > COa 4- H 




19 


O 4 


- D 


— > OD 




20 


O 4 


- DCO 


— > CO2 + D 




21 


OH 


+ D 


— > HDO 




22 


CO 


+ D 


— > DCO 


1930 


23 


HCO + D 


— > HDCO 




24 


H 2 CO + D 


— > H 2 DCO 


1799 


25 


H3CO + D 


— > CH3OD 




26 


D 4 


D 


— D 2 




27 


D 4 


OD 


— ► D2O 




28 


D 4 


- DCO 


— > D 2 CO 




29 


D 4 


- HDCO 


— > HD2CO 


1758 


30 


D 4 


- D 2 CO 


— > D3CO 


1713 


31 


D 4 


- H2DCO 


— > CH2DOD 




.32 


D 4 


- HD2CO 


— > CHD2OD 




.33 


D 4 


- D3CO 


— > CD3OD 





a See ICaselli et all J2002fl 

The surface reactions used in our model as well as any 
non-zero activation ene rgies j? a are given i n Table [3] More 
details can be found in lCaselli et alJ i2002T) . 



3 MASTER EQUATION FOR THE SYSTEM 

In the direct master equation method, the system corre- 
sponding to the grain surface is represented by a multitude 
of states that represent different possible populations of the 
surface species. In one state we might have particles of 
the species A, 1 particle of the species B, 2 particles of the 
species C, . . . , whereas a different state might consist of 1 
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particle of A, of B, 2 of C, ... . With each state, we asso- 
ciate a probability for the system to be in this state, and we 
develop equations for the rate at which these probabilities 
change (the master equation). Basic to the approach is the 
fact that the populations of some of the surface species are 
correlated. As a consequence, the method, at least formally, 
requires the consideration of the evolution of the system as 
a whole, as opposed to the rate equation approach, in which 
the evolution of the average population of any given species 
is followed separately. 

In the most general case, all the surface species are 
taken into consideration in each state of the system. Since 
the standard rate equations, however, give a sufficiently ac- 
curate description of the evolution of high-abundance species 
ijStantcheva et all2002h . the inclusion of these species in the 
master equation will impose an unnecessarily heavy load on 
the computing resources given the need to integrate many 
coupled equations simultaneously. Thus, a more practical 
approach involving fewer coupled equations is to include 
only the highly reactive and low-abundance surface species 
in the master equation, and to use rate-like equations to 
solve for the surface populations of the rest of the species. As 
a guide to determine which species should be included in the 
master equation, one can use the results obtained via rate 
equations and, if available, Monte Carlo simulations. More 
generally, these species are reactive atoms and radicals, es- 
pecially those that are precursors of major grain species. 
Here we include the 11 species H, O, D, OH, OD, HCO, 
DCO, CH 3 0, CH2DO, CHD2O, and CD 3 0. We enumerate 
these from 1 to 11 and refer to them as "probabilistic". The 
remainder of the species are referred to as normal and are 
treated via rate-like equations. 

Once we have determined the probabilistic species, we 
assign a probability P(ii, ■ . ■ , in) to any state {ii, . . . , in} 
that consists of ii H atoms, 12 O atoms, 13 D atoms, «4 OH 
molecules, . . . , in CD3O molecules, and solve the master 
equation, which consists of equations for the rate of change 
of the state probabilities P(ii, . . . , in). These equations are 
of the form 



,hi) = 



{X} 



fc acc (X)n(X) [P(...,ij 



+X)-! P (X) [fe + 1)-P(.-., h + 1, ■•■) - ijP(-,ij,-)] 

m 

{X,Y} (1) 

- ^2 kx,Y(ij)(ik)P(-,ij,-,ik,-) 

{X,Y} 



{X} 
{X} 



(i j + 2)(i j + l) 



P(...,i J +2, 



PC...,, 



of the species X, ik to that of the species Y, n(X) is 
the gas-phase abundance of X, and fc aC c(X), i^i p (X), and 
fcx,Y are the accretion and evaporation rates of X, and the 
rate coefficient for rea ction between X and Y ; respectively 
j jHaseeawa et al.ll992l) . In the units used here llCaselli et alJ 
1998), the rate coefficient is simply the sum of the the dif- 
fusion rates of X and Y over the entire grain. The first 
two terms on the right-hand side of eq. Q represent the 
changes in the probability with time due to accretion and 
evaporation processes. The subsequent terms account for the 
changes due to reactions between two probabilistic species, 
both non-identical and identical. It is important to note that 
there should be two more terms on the right-hand side of 
eq. Q, which are present in our calculations but are not 
given here for considerations of simplicity. The first of these 
terms takes into account the change of the probability due to 
reactions between a probabilistic and a normal species, such 
as H+CO — >HCO. The second accounts for (slow) reactions 
between two normal species, which have at least one proba- 
bilistic species as a product. In our current model there are 
no reactions of the latter type. In addition to these simpli- 
fications, we have also not indicated changes to more than 
two stochastic species in any term although there are re- 
actions in which the populations of three such species can 
change (e.g. O + HCO). 

The solution of eqs. (0 gives the probabilities for all the 
states, which, in turn, allows the calculation of the average 
abundance, or number of species per grain, of each proba- 
bilistic species, (iVx), as well as any necessary correlation 
terms, of the form (NxNy)- These entities, then, are substi- 
tuted into the equations used to solve for the abundances of 
the normal species. Although the latter equations are very 
similar to the widely used rate equations, they bear impor- 
tant differences. Equations (3J and © are the expressions 
for the average abundances of H2 and H2CO (both being 
normal species): 



dt 

d(NH 2 co) 
dt 



- <evip(H 2 ) X (N H2 ) 

+ 0.5fcH, H X (Nh(N h -1)), (2) 
+ fa,HCO x (A^hA^hco) 

- fc H ,H 2 CO X {Nil) X (iVHaCo). (3) 



It can be seen that, unlike regular rate equations, both 
correlated abundances and terms in which 1 is subtracted 
from an abundance appear. The rate-like equations are ob- 
tained by summation over single and correlated probabili- 



mation and depletion of major 


species 


Biham et alJl200ll: 


iGreen et al.ll200 ll IStantcheva et 


alJl200S 


)■ 



In eqs. (0, ij corresponds to the surface abundance 



To propagate the system forward in time, both the mas- 
ter equation and the rate-like equations must be solved si- 
multaneously. Since the probabilistic species have very low 
abundance, the probability for the system to be in states 
that contain high abundances of these species must be very 
small. It is reasonable, then, to neglect such states and 
limit our consideration to those that comprise only very few 
particles. For this purpose, we choose a set of parameters, 
{jVi,iV2,. • • , Nn} and iVtot, such that we neglect every state, 
{ii,i2,. . . ,in}, for which the following conditions hold true: 
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ij > Nj , for any j, or (4) 
^2 ij > iVtot • (5) 

3 

Because the first three particles in {ii,i2,- • • ,in} corre- 
spond to H, O, and D, we require that N 1} N2, and -/V3 be 
equal to at least 2 so that molecular hydrogen, oxygen and 
deuterium can be produced. Thus, the set of minimum val- 
ues for the upper limits Nj is {2,2,2,1,1,. . . ,1}. If we choose 
{Nj } to be equal to this set of minimum values, condition @ 
will eliminate all but TV = 3 3 x 2 8 = 6912 possible states 
of the system. This is still a very large number of coupled 
differential equations! A strong reduction in this number, 
however, can be obtained by applying condition which 
we have not done previously. Of course, the number of states 
will be dependent on the value of iVtot chosen. Clearly, any 
meaningful value of iVtot must be larger than or equal to 
each Nj and smaller than their sum, J^. Nj. If iVtot is taken 
to be 2, the number of states will be 70 for the set of limits 
{2,2,2,1,1,..., 1}. 

To understand how this number is determined, con- 
sider the following argument. Let the number of probabilistic 
particles be m. For the chosen Nj and iVtot, there will be 
four different groups of states: the first group, [0,0,. ..,0], 
consisting of only one state {0,0,. . . ,0}, the second group, 
[1,0,. . . ,0], consisting of m different states, the third group, 
[1,1,0,0,. ..,0], consisting of m ^,~ 1 ^ states, and the fourth 
group, [2,0,. . . ,0], consisting of only 3 states - {2,0,. . . ,0}, 
{0,2,0,. . . ,0}, and {0,0,2,0,. . . ,0}. In our model m=ll, so 
that the total number of states is 1 + 11 + 55 + 3 = 70. 

Although analytical formulae can be derived for the 
number of states determined by other specific values for m, 
{Ni,N%,. . . ,iV m } and iVtot, we choose to determine this num- 
ber by a computer subroutine. 



4 RESULTS 

Table0]shows a comparison between the surface populations 
(in monolayers per grain) of normal (stable, non-deuterated) 
species calculated via the direct master equation (ME) and 
the Monte Carlo (MC) approaches for the low, intermediate, 
and high density cases. To use the former method, all Nj 
were set to their minimum values except for N2 and N4, 
which correspond to O and OH. The calculations showed 
the latter two average abundances to be close to unity in 
some cases so that higher values for their upper limits had 
to be considered. 

Table 0] also shows the values of N2 and N4 for the 
three cases as well as the limit iVtot for the total number of 
probabilistic species and the total number of states M. One 
can see that for the high and intermediate density cases, the 
direct master equation method requires the solution of more 
than 816 coupled differential equations since the number 
M — 816 does not include the rate-like differential equations 
for the stable species. 

For the low and intermediate density cases, the agree- 
ment between the master equation and the Monte Carlo 
approaches is excellent, and the master equation approach 
outperforms the Monte Carlo calculation in terms of com- 
puter time, dramatically so for the low density limit. For 



the high density case, the agreement is not as excellent, but 
the results of the two methods deviate by at most 10%. This 
small deviation can be explained by the high abundance of O 
and OH on the surface, a condition that requires that higher 
values for N2 and N4, and consequently many more states, 
be considered in the master equation approach. In general, 
the agreement between the two formally exact approaches 
can be made as good as desired by considering sufficiently 
high values for the Nj and for 7V to t. This, however, will mean 
longer times for the master equation calculations to run, so 
one needs to find the right balance between the desired ac- 
curacy and the usage of computing resources. In general, it 
is a current weakness of the master equation approach that 
there is no obvious algorithm to determine the proper upper 
limit to the number of states to be considered. 

The fractionation ratios, /XD, are defined as the ra- 
tio between the abundances of the deuterated and the nor- 
mal isotopomers. The results for the fractionation ratios are 
given in Table |K| for the three cases discussed previously. 
The agreement between the two methods is reasonable but 
not perfect; here the major problem lies in the Monte Carlo 
approach, which is not accurate for small populations. In- 
deed, for some of the smaller fractionation ratios, no val- 
ues could be calculated with this method since populations 
less than unity are not allowed. As an example, consider 
the case of CHD2OH at intermediate density. From the ta- 
bles, one sees that the overall population of this molecule in 
the Monte Carlo approach is approximately 140 molecules. 
While the random (square-root) error is thus at the 10% 
level, the actual error is considerably greater and is proba- 
bly large enough to cover the 30-40% difference between the 
two methods of calculation. 

We have also used the direct master equation approach 
to model the chemistry in the temperature range 10-20 K 
for the high density case. Fig. shows the fractionation 
ratio / of doubly and multiply deuterated isotopomers in 
this temperature range, while a plot of the mole fraction 
of the normal isotopomers vs. temperature is given in Fig. 
[5] The parameters for the calculation (number of states, 
etc.) are the same as shown in Table 01 It can be seen that 
as the surface temperature increases from 10 to 20 K, the 
species (H2CO, CH3OH) made by hydrogenation involving 
reactions with activation energy decline sharply in abun- 
dance. This non-intuitive result derives from the fact that 
at higher temperatures, H atoms evaporate before they re- 
act. In addition, for most species, the fractionation ratios 
decline sharply with increasing temperature, especially for 
isotopomers that must be formed by reactions with activa- 
tion energy and involving D atoms, since these atoms tunnel 
more poorly than their H counterparts and are more likely 
to evaporate. 

Finally, the dependence of the fractionation ratios for 
doubly and multiply deuterated isotopomers on the ratio of 
atomic deuterium to hydrogen is shown in Fig. |3] for a 10 
K cloud at high density. The abscissa is actually the ratio 
of the accretion rate of atomic D to that of atomic H; this 
is equal to the ratio of the gas-phase abundances multiplied 
by a factor of 2 -0 ' 5 that derives from the relative speeds 
of D and H. Not surprisingly, the fractionation ratios in- 
crease as the D/H accretion ratio increases, with the slope 
proportional to the number of deuterium atoms in the iso- 
topomer. In fact, the results in Fig. [3] lie very close to the 
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Figure 1. Fractionation ratios / vs. T(j ust (K) at 10 4 yr for the high density case, with n(D)=0.3 cm -3 , computed with the direct master 
equation technique. The results lie close to a simple analytical limit near 10 K. 




Figure 2. Mole fractions of major surface species vs. T dust (K) at 10 4 yr for the high density case, computed with the direct master 
equation technique. 
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Table 4. Populations in mono-layers at 10 4 yr and 10 K. 



Species 


High Density 


Interm. 


Density 


Low Density 




MC 


ME 


MC 


ME 


MC 


ME 


CO 


5.0 


4.9 


0.00 


0.00 


0.0 


0.0 


H 2 


1.4 


1.2 


0.45 


0.45 


0.070 


0.069 


o 2 


2.7 


2.4 


0.080 


0.079 


0.0013 


0.0013 


co 2 


0.67 


0.71 


0.055 


0.055 


7.7(-4) 


7.2(-4) 


H 2 CO 


0.44 


0.44 


0.0 


0.0 


0.0 


0.0 


CH 3 OH 


0.088 


0.095 


0.42 


0.42 


0.046 


0.046 


Total abundance 


11.1 


10.5 


1.37 


1.36 


0.165 


0.164 


N 2 =N 




4 




4 




2 


A r 4=A f OH 




4 




4 




1 


Mot 




4 




4 




3 


N 




816 




816 




265 


CPU sec (Cray SV1) 


789 


313 


366 


245 


366 


14 



Table 5. Abundance ratios / at 10 4 yr and 10 K. 



Species High Density Interm. Density Low Density 





MC 


ME 


MC 


ME 


MC 


ME 


/HDCO 


0.33 


0.34 










P2CO 


0.026 


0.027 










iCH 3 OD 


0.20 


0.20 


0.19 


0.18 


0.18 


0.18 


yCH 2 DOH 


0.70 


0.71 


0.19 


0.19 


0.19 


0.19 


JCH2DOD 


0.14 


0.14 


0.034 


0.034 


0.034 


0.034 


/CHD2OH 


0.17 


0.18 


3.3(-4) 


2.1(-4) 


2.6(-4) 


1.4(-4) 


yCHD 2 OD 


0.034 


0.035 


6.4(-5) 


3.9(-5) 




2.5(-5) 


/CHD3OH 


0.014 


0.015 




5.6(-8) 




2.3(-8) 


/CHD3OD 


0.0025 


0.0029 




1.0(-8) 




4.3(-9) 


jHDO 


0.39 


0.38 


0.38 


0.38 


0.39 


0.39 


P2O 


0.038 


0.036 


0.036 


0.036 


0.036 


0.038 



very simple limit l|Tielensl 1198.4 iBrown fc Millarl Il98fll) in 
which all H and D atoms that accrete onto the surface of 
a grain eventually react with the CO reservoir despite the 
activation energy barriers. In this instance, if we let R be 
the ratio of the accretion rate of D to that of H, the fraction- 
ation ratios for the deuterated species are simple powers of 
R multiplied by a statistical factor, where the power is just 
the number of D atoms in the isotopomer. The statistical 
factor expresses the number of possible paths leading to the 
isotopomer, which, for example, is 3 for CH 2 DOH and 1 for 
CH3OD. This simple limit is independent of temperature, 
and so works progressively more poorly as the temperature 
is raised from 10 K, as can be seen from the temperature 
dependence in Fig. It also works more poorly for low and 
intermediate densities, where the CO reservoir does not ex- 
ist. 



5 DISCUSSION 

We have shown how a moderately complex system of sur- 
face chemical reactions can be solved by the formally exact 
direct master equation method. This method, or the equiva- 
lent Monte Carlo approach, is needed in the so-called accre- 
tion limit, which pertains to a situation in which the average 
abundance of important reactive species on a grain surface 



is less than unity. In the accretion limit, normal rate equa- 
tions are likely to be inaccurate, and a more exact method, 
which treats the discrete nature of the surface populations, 
is needed. Of the two exact methods used to date, the direct 
master equation approach offers the possibility of relative 
ease of implementation in combined gas-grain approaches 
since it involves the solutions of coupled differential equa- 
tions. 

In order to reach the level of complexity discussed here, 
we have introduced a new method to limit the number of 
states considered. This goal is achieved by limiting the to- 
tal number of reactive species considered on a grain. In the 
current system, which involves 11 reactive species, we have 
limited consideration to a total number of species in the 
range 3-4. Of course, this approach only works successfully 
if the probability for even this number of total species is 
small, which is what defines the accretion limit. 

If one eliminates D and the five deuterium-containing 
radicals from the 11 reactive species considered here, and 
replaces them with other reactive species leading to major 
grain constituents, it is entirely possible that the master 
equation treatment described here can be extended to study 
combined gas-grain models of considerable complexity, in 
which gas-phase abundances change with time. Problems 
to be solved before this goal becomes a reality include the 
need to develop a criterion for the smallest number of states 
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Figure 3. Fractionation ratios / for high density conditions at 10 K vs. acc(D)/acc(H, the ratio of the accretion rate of atomic D to 
atomic H. 



that should be considered, and the related need to determine 
when and for which species normal rate equations can be 
utilized. These problems are certain to depend on physical 
conditions such as grain size and temperature, as well as 
the rates of diffusion chosen. It would be ideal if a program 
could be developed to decide these issues as the integration 
with time proceeds. 



ACKNOWLEDGMENTS 

We thank the National Science Foundation (US) for support 
of the Ohio State program in astrochemistry and the Ohio 
Supercomputer Center for time on their Cray SV1 machine. 
Special thanks go to V. I. Shematovich for carefully reading 
a version of the manuscript. 



Of the deuterated species discussed here, HDCO, 
D 2 CO, CH3OD, CH2DOH, CHD2OH, and HDO have al- 
ready been detected in interstellar sources. A low-mass pro- 
tostellar source - IRAS 16293-2422 - contains both deuter- 
ated isotopomers of formaldehyde, both singly deuterated 
isotopomers of me thanol, and CHD2OH in the gas phase 
JParise et alJl2002j) . It is likely that at least some of the 
fractionation occurs on interstell ar grains and is follo wed 
by evapo ration back into the gas dCharnlev et al]ll997f) . As 
shown bv lParise et al.l (|2002); howe ver, the granu l ar deu ter- 
ation model presented her e and in | Casell i et all i2002T) . as 
well as the simpler one in ICharnle^t^dTi^gTl) . does not 
account quantitatively for all of the observations with any 
single value of gas-phase D chosen. Clearly, a more complex 
gas-grain model is needed, in which time-dependent fraction- 
ation occurs both in the gas and on dust-particle surfaces. 
Our current success with the direct master equation method 
may well allow us to consider such a complex gas-grain sys- 
tem in the future. 
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